#ifndef _PP_FORM_
#define _PP_FORM_
#include "Function.h"
#include "spline.h"
#include <vector>
#include <eigen3/Eigen/Sparse>
#include <eigen3/Eigen/Dense>
#include <iostream>
using namespace std;

class Polynomial : public Function
{
private:
    vector<double> co;
    double x_0;
    int n;
public:
    Polynomial()=default;
    Polynomial(vector<double> _co, double _base): co(_co), x_0(_base) {}
    double operator()(double _x)
    {
        double tmp = co[0];
        n = co.size();
        for (int i = 1; i < n; i++)
        {
            double delta_tmp = co[i];

            for (int j = 0; j < i; j++)
            {
                delta_tmp = delta_tmp*(_x - x_0);
            }

            tmp = tmp + delta_tmp;
        }
        return tmp;
    }
};



class PP_Form : public spline
{
private:
    vector<double> point;
    Function &f;
    vector<Polynomial> poly; // 存储插值得到的n-1个分段多项式
    vector<vector<double> > Co;
    int n, choice, k; 
    Eigen::VectorXd co_m; // co_m：解方程最终得到的系数(一阶导数m)
public:
    // f:待插值函数；point：插值点t_1,...t_n；
    // choice：插值边界条件（1为CCS, 2为CSS, 3为NCS）；
    // k：插值样条基函数阶数。
    PP_Form(Function &_f, vector<double> _point, int _choice, int _k): f(_f), point(_point), choice(_choice), k(_k)
    {
        if ((choice != 1) && (choice != 2) && (choice != 3))
	    {
		    cerr<< "error: the choice is not available." <<endl;
		    exit(-1);
	    }
        if ((k != 1) && (k != 3))
	    {
		    cerr<< "error: the k is not available." <<endl;
		    exit(-1);
	    }

        n = point.size();
    }

    void cal()
    {
        if (k == 3)
        {
            Eigen::SparseMatrix<double> A(n, n);
            vector<Eigen::Triplet<double> > tripletlist;
            Eigen::MatrixXd y(n, 1); 
            y = Eigen::MatrixXd::Zero(n, 1); // 解方程Am = y

            for (int i = 1; i < n - 1; i++)
		    {
                tripletlist.push_back(Eigen::Triplet<double>(i, i-1, (point[i+1]-point[i])/(point[i+1]-point[i-1])));
                tripletlist.push_back(Eigen::Triplet<double>(i, i, 2.0));
                tripletlist.push_back(Eigen::Triplet<double>(i, i+1, (point[i]-point[i-1])/(point[i+1]-point[i-1])));
                y(i, 0) = 3.0*(point[i]-point[i-1])/(point[i+1]-point[i-1])*(f(point[i+1])-f(point[i]))/(point[i+1]-point[i]) + 3*(point[i+1]-point[i])/(point[i+1]-point[i-1])*(f(point[i])-f(point[i-1]))/(point[i]-point[i-1]);
		    }

            if (choice ==  1)
            {
                tripletlist.push_back(Eigen::Triplet<double>(0, 0, 1.0));
                tripletlist.push_back(Eigen::Triplet<double>(n-1, n-1, 1.0));
                y(0, 0) = f.diff(point[0]);
                y(n-1, 0) = f.diff(point[n-1]);
            }
            else
            {
                tripletlist.push_back(Eigen::Triplet<double>(0, 0, 4.0));
                tripletlist.push_back(Eigen::Triplet<double>(0, 1, 2.0));
                tripletlist.push_back(Eigen::Triplet<double>(n-1, n-2, 2.0));
                tripletlist.push_back(Eigen::Triplet<double>(n-1, n-1, 4.0));
                if (choice ==  2)
                {
                    y(0, 0) = 6.0*(f(point[1])-f(point[0]))/(point[1]-point[0]) - f.diff2(point[0])*(point[1]-point[0]);
                    y(n-1, 0) = 6.0*(f(point[n-1])-f(point[n-2]))/(point[n-1]-point[n-2]) + f.diff2(point[n-1])*(point[n-1]-point[n-2]);
                }
                if (choice ==  3)
                {
                    y(0, 0) = 6.0*(f(point[1])-f(point[0]))/(point[1]-point[0]);
                    y(n-1, 0) = 6.0*(f(point[n-1])-f(point[n-2]))/(point[n-1]-point[n-2]);
                }
            }

            A.setFromTriplets(tripletlist.begin(), tripletlist.end());
            A.makeCompressed();
            Eigen::SparseLU<Eigen::SparseMatrix<double> > cal_co;
            cal_co.compute(A);
            co_m = cal_co.solve(y);

            for (int i = 1; i < n; i++)
		    {
                vector<double> pp;
                double K; 
                K = (f(point[i])-f(point[i-1]))/(point[i]-point[i-1]);
                pp.push_back(f(point[i-1]));
                pp.push_back(co_m(i-1));
                pp.push_back((3.0*K-2.0*co_m(i-1)-co_m(i))/(point[i]-point[i-1]));
                pp.push_back((co_m(i-1)+co_m(i)-2.0*K)/pow(point[i]-point[i-1],2));
                Co.push_back(pp);

                Polynomial p(pp, point[i-1]);
                poly.push_back(p); 
		    }
        }
        if (k == 1)
        {
            for (int i = 1; i < n; i++)
		    {
                vector<double> pp;
                double K; 
                K = (f(point[i])-f(point[i-1]))/(point[i]-point[i-1]);
                pp.push_back(f(point[i])-K*point[i]);
                pp.push_back(K);
                Co.push_back(pp);
                
                Polynomial p(pp, 0);
                poly.push_back(p); 
		    }
        }
    }

    double operator()(double _x)
    {
        if ((_x < point[0]) || (_x > point[n-1]))
	    {
		    cerr<< "error: ut of range." <<endl;
		    exit(-1);
	    }
        if (_x == point[0])
         {   return f(point[0]);}

        int position = 1;
        while (_x > point[position])
        {
            position = position + 1;
        }

        if (_x == point[position])
            return f(point[position]);
        else
            return poly[position-1](_x);
    }

};

#endif